close all
clear
clc
%%%%%%%%%%%%%%%%%%%%%
%     Parameter     %
%%%%%%%%%%%%%%%%%%%%%
a=0.2;                      %lowpass filter not used
slideWindow=5;             %delta function forwarding and tailing window
cutoff=15;                  %peak selection cutoff critias
excludePercentage=0;      %fitting excluding percentage
fitWindow=15;
cutoffForce=9;
stepCutoffMagnet=9;
internalTTest=1;            %Single T test
internalTTestAlpha=0.01;
position_overstretching=12.41;
aa_length=0.5;    % the pesistent length of a amino acid
filelist=importdata('G:\catenin\wt\covhis_sur\filelist_vh.txt');
only_unfolding=1;% only process the data in the unfolding part of spectrum
with_smoothing=0;%whether to smooth rawdata
drift_correction=0; %whether to correct linear drift in long time measurements

%change to dump if the analysis is repeated to avoid double counting
summaryoutputfile=('G:\catenin\wt\covhis_sur\goodunfolding_vh.txt');
good_threshold=50;
fid1=fopen(summaryoutputfile,'a');


for currentfile = filelist'
    
    name=char(currentfile);
    
    
    %% Data processing
    steps_data=[];
    filename=sprintf('%s%s',name,'.txt')
    %c=65/(exp(-(13.5-position_overstretching)/0.36)+0.48*exp(-(13.5-position_overstretching)/1.12));
    c=280;
    datastruct=importdata(filename);
    dataraw=datastruct.data;
    egxT=dataraw(:,1)';
    egMagneticT=dataraw(:,2);
    egRawDynT=dataraw(:,3);
    %plot(egxT); % Check equal spacing
    
    %Mingxi data: A. linear drift remover
    %             B. force altering regions
    %             C. points for fitting selection
    
    if drift_correction
        linearFitP=polyfit(egxT(egMagneticT==min(egMagneticT))',egRawDynT(egMagneticT==min(egMagneticT)),1);
        egfittedDriftT=polyval(linearFitP, egxT)';
        disp('Determining the drift\n');
        linearFitP % display the drift
        egRawDynTShift=egRawDynT-egfittedDriftT;
    else
        egRawDynTShift=egRawDynT;
    end
    %figure;
    %plot(egxT,egRawDynT);
    %hold on; plot(egxT,egfittedDriftT,'r','LineWidth',3);
    %figure;
    %plot(egxT,egRawDynT-egfittedDriftT,'b');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %cutoffForce=min(egMagneticT);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %process if only unfolding signal is processed
    tsig=egMagneticT>cutoffForce;
    dsig=diff([0;tsig;0]);
    startIndex=find(dsig > 0);
    endIndex=find(dsig < 0)-2;
    
    
    hFig = figure();
    set(hFig, 'Visible', 'off')
    
    
    for forceAlteringRegionIndex=1:length(startIndex)
        is_good=0;% whether this tether has enough unfolding steps
        single_steps_data=[];
        single_cumulative_contour=0; % total unfolding size in a single force altering region
        region_end_index=endIndex(forceAlteringRegionIndex);
        region_start_index=startIndex(forceAlteringRegionIndex);
        if only_unfolding    %with this option only the signals in the unfolding region of the trajectory is analyzed
            region_Magnetic=egMagneticT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex)); % find the magnetic positions in the region
            region_magsig=diff([0;region_Magnetic;0]);
            if ~isempty(find(region_magsig<0,1,'first'))
                region_end_index=region_start_index+find(region_magsig<0,1,'first')-1;
            end
        end
        %region_start_index
        %region_end_index
        egx=egxT(region_start_index:region_end_index);
        egRawDyn=egRawDynTShift(region_start_index:region_end_index);
        %scaled magnetic Location
        egMagnetic=egMagneticT(region_start_index:region_end_index);
        
        %egTruDyn=(egMagnetic-min(egMagnetic))*(max(egRawDyn)-min(egRawDyn))/(max(egMagnetic)-min(egMagnetic));
        %figure;
        %plot(egxT,egRawDynTShift);
        %hold on;plot(egx,egRawDyn,'r');
        
        
        %% Assign
        x=egx;
        RawDyn=egRawDyn;
        %TruDyn=egTruDyn;
        %% Filtering By lowpass filter (Freqency**)
        % Without filter seems better for loading rate data (or change parameter)
        FilDynA=filter(a,[1 a-1],RawDyn);
        FilDynB=smooth(RawDyn,slideWindow,'sgolay'); % SG smooth
        %
        
        subplot(3,1,1);
        plot(x,RawDyn,'b',x,FilDynA,'r');
        hold on; plot(x,FilDynB,'g');
        
        
        %% auxiliary function (delta function)
        
        disp([num2str(slideWindow), ' points delta function ...']);
        if with_smoothing
            customDyn=FilDynB;
        else
            customDyn=RawDyn;
        end
        deltaF=zeros(1,numel(customDyn));
        
        %starting
        for i=1:slideWindow
            diff_value=i-1;
            forwarding=mean(customDyn(i:i+diff_value));
            tailing=mean(customDyn(i-diff_value:i));
            deltaF(i)=tailing-forwarding;
        end
        
        for i=slideWindow+1:numel(customDyn)-slideWindow
            %disp(i);
            forwarding=mean(customDyn(i:i+slideWindow));
            tailing=mean(customDyn(i-slideWindow:i));
            deltaF(i)=tailing-forwarding;
        end
        
        %ending
        for i=numel(customDyn)-slideWindow+1:numel(customDyn)
            diff_value=numel(customDyn)-i;
            forwarding=mean(customDyn(i:i+diff_value));
            tailing=mean(customDyn(i-diff_value:i));
            deltaF(i)=tailing-forwarding;
        end
        
        
        %% peak location determination
        
        disp([num2str(cutoff), ' cutoff peak locating ...']);
        
        deltaFCut=deltaF.*(deltaF>=cutoff|deltaF<=-cutoff);
        deltaFCutU=deltaF.*(deltaF>=cutoff);  %positive
        deltaFCutD=deltaF.*(deltaF<=-cutoff); %negative
        
        
        [Upeaks,Ulocs]=findpeaks(deltaFCutU,'minpeakdistance',slideWindow);
        [Dpeaks,Dlocs]=findpeaks(-deltaFCutD,'minpeakdistance',slideWindow);
        
        subplot(3,1,3);
        plot(x,deltaF,'k-');
        hold on;plot(x,deltaFCut,'g');
        plot(Ulocs*(max(x)-min(x))/(numel(x)-1)+min(x),Upeaks,'ro');
        plot(Dlocs*(max(x)-min(x))/(numel(x)-1)+min(x),-Dpeaks,'mo');
        
        %%%%%%%%%%%%%%%%%%%%%%
        %        TRAIL       %
        %%%%%%%%%%%%%%%%%%%%%%
        %% T Test filter remove false steps
        fitDyn=RawDyn;
        locs=sort([1,Ulocs, Dlocs, numel(fitDyn)]);
        
        % Multiple rounds of test until fix
        if internalTTest == 1
            compareDyn=fitDyn;
            compareDyn_size=size(compareDyn);
            templocs=zeros(1,numel(locs));
            previouslocs=0;
            count=1;
            
            while numel(previouslocs) ~= numel(locs)
                disp([num2str(internalTTestAlpha*100), '% Internal Student T Test: ', num2str(count)]);
                previouslocs = locs;
                for i=2:numel(locs)-1
                    %interval selection
                    startI=locs(i)-fitWindow;
                    endI=locs(i)+fitWindow;
                    if startI < 1
                        startI=1;
                    end
                    if endI > compareDyn_size
                        endI=compareDyn_size;
                    end
                    
                    jumpI=locs(i);
                    startCutI=startI+floor((jumpI-startI)*excludePercentage);
                    endCutI=endI-floor((endI-jumpI)*excludePercentage);
                    groupA=compareDyn(startCutI:jumpI-1);
                    groupB=compareDyn(jumpI+1:endCutI);
                    
                    [H,PI]=ttest2(groupA,groupB,internalTTestAlpha,[],'unequal');
                    %disp(PI);
                    
                    if H==1
                        templocs(i)=locs(i);
                    else
                        locs(i)=startI;
                    end
                end
                
                locs=[1,templocs(templocs~=0), numel(fitDyn)];
                templocs=zeros(1,numel(locs));
                count=count+1;
            end
        end
        
        %% least square fit to heaviside function (per interval)
        %Change to fit locally to adjacentpoitns
        fittedDyn=zeros(numel(fitDyn),1);
        
        for i=2:numel(locs)-1
            %tempfix for some locs<0 error
            if locs(i)>fitWindow && locs(i)<(length(x)-fitWindow)
                %interval selection
                startI=locs(i)-fitWindow;
                endI=locs(i)+fitWindow;
                jumpI=locs(i);
                startCutI=startI+floor((jumpI-startI)*excludePercentage);
                endCutI=endI-floor((endI-jumpI)*excludePercentage);
                currentX=x(startCutI:endCutI);
                currentY=fitDyn(startCutI:endCutI);
                
                %figure;plot(x(startI:endI),fitDyn(startI:endI),'go',currentX,currentY,'ro')
                %guess starting points
                ga=mean(fitDyn(jumpI:endCutI))-mean(fitDyn(startCutI:jumpI));
                gb=x(jumpI);
                if ga>=0
                    gc=mean(fitDyn(startCutI:jumpI));
                else
                    gc=mean(fitDyn(jumpI:endCutI));
                end
                gd=polyfit(x(startCutI:jumpI)',fitDyn(startCutI:jumpI),1);
                %fit (added a linear term to approximate the loading rate
                fequation='a.*heaviside(x-b)+c+d.*(x-b)';
                ft = fittype(fequation);
                fo = fitoptions('Method','NonlinearLeastSquares','StartPoint',[ga gb gc gd(1)],'Lower',[-inf min(x) -inf -inf],'Upper',[inf max(x) inf inf]);
                currentFit = fit(currentX',currentY,ft,fo);
                
                %currentF=currentFit(currentX);
                %hold on;plot(currentX,currentF,'k','LineWidth',3);
                %currentFit
                
                %update result
                locs(i)=round((currentFit.b-min(x))/(max(x)-min(x))*(numel(x)-1));
                
                leftValue=currentFit.c;
                rightValue=currentFit.a+currentFit.c;
                stepValue=0.5*currentFit.a+currentFit.c;
                if locs(i)>fitWindow
                    fittedDyn((locs(i)-fitWindow):locs(i)-1)=leftValue;
                else
                    fittedDyn(1:locs(i)-1)=leftValue;
                end
                fittedDyn((locs(i)+1):locs(i)+fitWindow)=rightValue;
                fittedDyn(locs(i))=stepValue;
                %get the parameters for the magnet and calculated the force
                stepsize=currentFit.a;
                locmagnet=egMagnetic(locs(i));
                force_then=c*(exp(-(13.5-locmagnet)/0.36)+0.48*exp(-(13.5-locmagnet)/1.12));
                full_length=stepsize/zz_wlc(force_then*aa_length/4);
                %record data for forces > preset values to reduce false positives
                if locmagnet>stepCutoffMagnet
                    steps_data=[steps_data;[forceAlteringRegionIndex,locmagnet,force_then,stepsize,full_length]];
                    single_cumulative_contour=single_cumulative_contour+full_length;
                    single_steps_data=[single_steps_data;[forceAlteringRegionIndex,locmagnet,force_then,stepsize,full_length]];
                end
            end
            
        end
        
        
        if single_cumulative_contour>=good_threshold
            is_good=1;
        end
        
        %write to global file if the tether is good
        if is_good
            
            for row= 1:size(single_steps_data,1)
                fprintf(fid1,'%s\t%d\t%.3f\t%.3f\t%.3f\t%.3f\n',name,single_steps_data(row,:)');
            end
        end
       
        % %update last step & first value
        % rightValue=currentFit.c+currentFit.a;
        % fittedDyn(locs(i)+1:numel(fittedDyn))=rightValue;
        % fittedDyn(1)=fittedDyn(2);
        
        % %update locs value
        % for i=2:numel(locs)-1
        %     fittedDyn(locs(i))=.5*(fittedDyn(locs(i)-1)+fittedDyn(locs(i)+1));
        % end
        
        subplot(3,1,2);
        plot(x,fitDyn,'ro');
        hold on; plot(x,fittedDyn,'g','LineWidth',3);
    end
    figure_filename = sprintf('%s%s',name,'.png');
    
    print(gcf, '-dpng','-opengl', figure_filename, '-r640');
    %open write file
    if only_unfolding
        filename = sprintf('%s%s',name,'-f-unfoldingstep.dat');
    else
        filename = sprintf('%s%s',name,'-f-step.dat');
    end
    fid = fopen(filename,'w');
    if(fid == 0)
        disp('file cannot be opened');
    end
    fprintf(fid,'%s\n',name);
    fprintf(fid,'rampnum\tmagnet(mm)\tforce (pN)\tstepsize(nm)\tconter length(nm)\n');
    fprintf(fid,'%d\t%.3f\t%.3f\t%.3f\t%.3f\n',steps_data');
    fclose(fid);
    
    
end

fclose(fid1);